source("regression_functions.R")
fig.path <- "results/"



### define outcomes and control variables
yvars<- paste0("exogamy_l",c(1:11,15:16)) ## 
yvars_std <- paste0(yvars,"_std")
controls.geo <- c("ruggedness_nunn_puga" ,  "elevation"  ,  "malaria_suit_max",
                  "TSI_CRU_mean_1901_1920" , "temperature_fao" , "precipitation_fao",  "coast_log"  ,
                  "LONGNUM", "agric_suit", "LATNUM")
controls.hist <- c("explorers_log", "cities_log", "capital_log","dist.prot_log" , "dist.print_log","popc_mean_1720_1890_log")
controls.ind <- c("age","age.sq","age.m","age.m.sq")
base.x <-  c(controls.ind, controls.geo,controls.hist)

### define analysis countries
sample_countries_hance <- unique(dhs$iso3c[dhs$sample_country_hance==1])


#### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### 
##### local ethnic diversity: DHS-based ELF scores at level of survey cluster
##### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### 


med.elf <- median(dhs$elf_cluster_side[!is.na(dhs$pubspc.23_poly) & dhs$LONGNUM!=0  & 
                                         dhs$iso3c%in%sample_countries_hance & !is.na(dhs$exogamy_l11)],na.rm=T)

dhs$elf_cluster_high <- ifelse(dhs$elf_cluster_side>med.elf,1,0)

m.list.ia <- lapply(c(1:13), function(i){
  felm(as.formula(RegFor( y = yvars_std[i] , x = c("hance_crops5_sum_15km_std*elf_cluster_high","pubspc.23_poly_std*elf_cluster_high", base.x) ,
                          FE = "country_survey_round" , 
                          IV="0", clust = "loc.id_pubs_poly")),
       data=dhs,subset=LONGNUM!=0 & iso3c%in%sample_countries_hance)
})


plot.df <- do.call(rbind,lapply(m.list.ia,function(m){
  rbind(
    deltaMethod(m,"hance_crops5_sum_15km_std"),
    deltaMethod(m,"hance_crops5_sum_15km_std + `hance_crops5_sum_15km_std:elf_cluster_high`"),
    deltaMethod(m,"pubspc.23_poly_std"),
    deltaMethod(m,"pubspc.23_poly_std + `elf_cluster_high:pubspc.23_poly_std`")
  )
}))


ticks <- rev(seq(1,19.5,1.5))

plot.df$dv <- rep(c(paste0("Exogamy L",c(1:10)),paste0("Exogamy L","11-14"),paste0("Exogamy L",c(15:16))),each=4)
plot.df$what <- rep(c("Crops (ELF Low)","Crops (ELF High)","Pub. (ELF Low)","Pub. (ELF High)"),13)

ord. <- c(ticks + 0.4,ticks + 0.2,ticks-0.2,ticks-0.4)
plot.df$order <- ord.[order(ord.,decreasing=T)]

plot.df$what <- factor(plot.df$what,levels=c("Crops (ELF Low)","Crops (ELF High)","Pub. (ELF Low)","Pub. (ELF High)"))


col1 <- wes_palette("Darjeeling1",5,"discrete")[2]
col2 <- wes_palette("Darjeeling1",5,"discrete")[3]
col3 <- wes_palette("Darjeeling1",5,"discrete")[5]
col4 <- wes_palette("Darjeeling1",15,"continuous")[2] 

names(plot.df)[c(1:4)] <- c("coef","se","CI_lower", "CI_upper")

# do the plot
p1 <- ggplot(plot.df)
p1 <- p1 + geom_point(size=2.25,aes(x=coef,y=order,color=what,shape=what,fill=what)) + 
  geom_errorbarh(aes(x=coef,y=order, xmin=CI_lower, xmax=CI_upper,color=what), size=1, height=0.0) +
  geom_vline(xintercept=0,linetype="dotted", size=0.6) +
  scale_y_continuous(breaks=ticks,minor_breaks = NULL, 
                     labels=c(paste0("Exogamy L",c(1:10)),"Exogamy L11-14",paste0("Exogamy L",c(15:16)))) + 
  #scale_x_continuous(minor_breaks = seq(0 , 0.55, 0.05), breaks = seq(0, 0.55, 0.1)) +
  labs(x = "Marginal Effects & 95% CIs", y = "Exogamy at Ethnologue Levels",
       title="Interactive Effects?",
       subtitle="Treatments interacted with local-level diversity"
       #,caption="Cash crop value per sqkm in 1960 USD (Poly.)\nCash Crops (SAR)\nCash Crop Suitability (Poly)\nCash Crops (Instrumented with Suitability)"
  ) +
  theme_minimal(base_size=14) +  
  theme(axis.text.y = element_text(size=10),legend.position = "bottom") +
  scale_color_manual(values=c("Crops (ELF Low)" = col1, "Crops (ELF High)" = col2, "Pub. (ELF Low)" =col3, "Pub. (ELF High)" =col4),name="Marginal Effect") + 
  scale_fill_manual(values=c("Crops (ELF Low)" = col1, "Crops (ELF High)" = col2, "Pub. (ELF Low)" =col3, "Pub. (ELF High)" =col4),name="Marginal Effect") +
  scale_shape_manual(values=c("Crops (ELF Low)" = 24, "Crops (ELF High)" = 22, "Pub. (ELF Low)" = 23, "Pub. (ELF High)" =25),name="Marginal Effect")
p1 <- p1 + guides(color=guide_legend(ncol=2,nrow=2,byrow=TRUE),
                fill=guide_legend(ncol=2,nrow=2,byrow=TRUE),
                shape=guide_legend(ncol=2,nrow=2,byrow=TRUE))
p1




#### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### 
##### local ethnic diversity: WLMS-based ELF scores at level of 05 degree grid.cells
##### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### ##### 


med.elf <- median(dhs$terr_frac_cell_alt[!is.na(dhs$pubspc.23_poly) & dhs$LONGNUM!=0  & 
                                         dhs$iso3c%in%sample_countries_hance & !is.na(dhs$exogamy_l11)],na.rm=T)

dhs$terr_frac_cell_high <- ifelse(dhs$terr_frac_cell_alt>med.elf,1,0)

m.list.ia <- lapply(c(1:13), function(i){
  felm(as.formula(RegFor( y = yvars_std[i] , x = c("hance_crops5_sum_15km_std*terr_frac_cell_high",
                                                   "pubspc.23_poly_std*terr_frac_cell_high", base.x) ,
                          FE = "country_survey_round" , 
                          IV="0", clust = "loc.id_pubs_poly")),
       data=dhs,subset=LONGNUM!=0 & iso3c%in%sample_countries_hance)
})


plot.df <- do.call(rbind,lapply(m.list.ia,function(m){
  rbind(
    deltaMethod(m,"hance_crops5_sum_15km_std"),
    deltaMethod(m,"hance_crops5_sum_15km_std + `hance_crops5_sum_15km_std:terr_frac_cell_high`"),
    deltaMethod(m,"pubspc.23_poly_std"),
    deltaMethod(m,"pubspc.23_poly_std + `terr_frac_cell_high:pubspc.23_poly_std`")
  )
}))


ticks <- rev(seq(1,19.5,1.5))

plot.df$dv <- rep(c(paste0("Exogamy L",c(1:10)),paste0("Exogamy L","11-14"),paste0("Exogamy L",c(15:16))),each=4)
plot.df$what <- rep(c("Crops (ELF Low)","Crops (ELF High)","Pub. (ELF Low)","Pub. (ELF High)"),13)

ord. <- c(ticks + 0.4,ticks + 0.2,ticks-0.2,ticks-0.4)
plot.df$order <- ord.[order(ord.,decreasing=T)]

plot.df$what <- factor(plot.df$what,levels=c("Crops (ELF Low)","Crops (ELF High)","Pub. (ELF Low)","Pub. (ELF High)"))


col1 <- wes_palette("Darjeeling1",5,"discrete")[2]
col2 <- wes_palette("Darjeeling1",5,"discrete")[3]
col3 <- wes_palette("Darjeeling1",5,"discrete")[5]
col4 <- wes_palette("Darjeeling1",15,"continuous")[2] 

names(plot.df)[c(1:4)] <- c("coef","se","CI_lower", "CI_upper")

# do the plot
p2 <- ggplot(plot.df)
p2 <- p2 + geom_point(size=2.25,aes(x=coef,y=order,color=what,shape=what,fill=what)) + 
  geom_errorbarh(aes(x=coef,y=order, xmin=CI_lower, xmax=CI_upper,color=what), size=1, height=0.0) +
  geom_vline(xintercept=0,linetype="dotted", size=0.6) +
  scale_y_continuous(breaks=ticks,minor_breaks = NULL, 
                     labels=c(paste0("Exogamy L",c(1:10)),"Exogamy L11-14",paste0("Exogamy L",c(15:16)))) + 
  #scale_x_continuous(minor_breaks = seq(0 , 0.55, 0.05), breaks = seq(0, 0.55, 0.1)) +
  labs(x = "Marginal Effects & 95% CIs", y = "Exogamy at Ethnologue Levels",
       title="Interactive Effects?",
       subtitle="Treatments interacted with local-level diversity"
       #,caption="Cash crop value per sqkm in 1960 USD (Poly.)\nCash Crops (SAR)\nCash Crop Suitability (Poly)\nCash Crops (Instrumented with Suitability)"
  ) +
  theme_minimal(base_size=14) +  
  theme(axis.text.y = element_text(size=10),legend.position = "bottom") +
  scale_color_manual(values=c("Crops (ELF Low)" = col1, "Crops (ELF High)" = col2, "Pub. (ELF Low)" =col3, "Pub. (ELF High)" =col4),name="Marginal Effect") + 
  scale_fill_manual(values=c("Crops (ELF Low)" = col1, "Crops (ELF High)" = col2, "Pub. (ELF Low)" =col3, "Pub. (ELF High)" =col4),name="Marginal Effect") +
  scale_shape_manual(values=c("Crops (ELF Low)" = 24, "Crops (ELF High)" = 22, "Pub. (ELF Low)" = 23, "Pub. (ELF High)" =25),name="Marginal Effect")
p2 <- p2 + guides(color=guide_legend(ncol=2,nrow=2,byrow=TRUE),
                  fill=guide_legend(ncol=2,nrow=2,byrow=TRUE),
                  shape=guide_legend(ncol=2,nrow=2,byrow=TRUE))
p2

# combine with patchwork
p.out <- p1 | p2
#save
p.out + ggsave(paste0(fig.path,"dhs_figA35.pdf"),width=14,height=7.5)
